3PI algorithm for spiral CT

ABSTRACT

Methods and systems for reconstructing images of moving objects being spirally scanned with two dimensional detectors with a 3PI algorithm. The moving objects can be scanned at a rate of up to approximately three times slower than those of pre-existing systems. In a preferred embodiment, the invention allows for a patient on a table moving through a spiral scanner to be slowed down by a factor of up to three, and still use the same size detector array as those in existing spiral scanning systems.

[0001] This invention claims the benefit of priority to U.S. ProvisionalApplication Serial No. 60/430,802 filed, Dec. 4, 2002, and is aContinuation-In-Part of U.S. patent application Ser. No. 10/389,534filed Mar. 14, 2003 which is a Continuation-In-Part of Ser. No.10/389,090 filed March Continuation-In-Part of Ser. No. 10/143,160 filedMay 10, 2002 now U.S. Pat. No. 6,574,299, which claims the benefit ofpriority to U.S. Provisional Application 60/312,827 filed Aug. 16, 2001.

FIELD OF INVENTION

[0002] This invention relates to computer tomography, and in particularto processes, methods and systems for reconstructing three-dimensionalimages from the data obtained by spiral scans using 3PI algorithm.

BACKGROUND AND PRIOR ART

[0003] Over the last thirty years, computer tomography (CT) has gonefrom image reconstruction based on scanning in a slice-by-slice processto spiral scanning. From the 1970s to 1980s the slice-by-slice scanningwas used. In this mode the incremental motions of the patient on thetable through the gantry and the gantry rotations were performed oneafter another. Since the patient was stationary during the gantryrotations, the trajectory of the x-ray source around the patient wascircular. Pre-selected slices through the patient have beenreconstructed using the data obtained by such circular scans. From themid 1980s to present day, spiral type scanning has become the preferredprocess for data collection in CT. Under spiral scanning a table withthe patient continuously moves through the gantry that is continuouslyrotating about the table. At first, spiral scanning has usedone-dimensional detectors, which receive data in one dimension (a singlerow of detectors). Later, two-dimensional detectors, where multiple rows(two or more rows) of detectors sit next to one another, have beenintroduced. In CT there have been significant problems for imagereconstruction especially for two-dimensional detectors. In what followsthe data provided by the two-dimensional detectors will be referred toas cone-beam (CB) data or CB projections.

[0004] For three-dimensional (also known as volumetric) imagereconstruction from the data provided by a spiral scan withtwo-dimensional detectors, there are two known groups of algorithms:Exact algorithms and Approximate algorithms, that each have knownproblems. Under ideal circumstances, exact algorithms can provide areplication of an exact image. Thus, one should expect that exactalgorithms would produce images of good quality even under non-ideal(that is, realistic) circumstances. However, exact algorithms can beknown to take many hours to provide an image reconstruction, and cantake up great amounts of computer power when being used. Thesealgorithms can require keeping considerable amounts of cone beamprojections in memory. Additionally, some exact algorithms can requirelarge detector arrays to be operable and can have limits on the size ofthe patient being scanned.

[0005] Approximate algorithms possess a filtered back projection (FBP)structure, so they can produce an image very efficiently and using lesscomputing power than Exact algorithms. However, even under the idealcircumstances they produce an approximate image that may be similar tobut still different from the exact image. In particular, Approximatealgorithms can create artifacts, which are false features in an image.Under certain circumstances these artifacts could be quite severe.

[0006] The first group of algorithms that are theoretically exact andhave the shift-invariant FBP structure was disclosed in U.S. patentapplication Ser. No. 10/143,160 filed May 10, 2002, now U.S. Pat. No.6,574,299, which claims the benefit of priority to U.S. ProvisionalApplication 60/312,827 filed Aug. 16, 2001. A shortcoming of thesealgorithms is that they operate in the 1PI mode. This means that if thedetector array is large in the axial direction, then one has totranslate the patient through the gantry very quickly in order to useall of the detector. However, rapid motion of the patient is veryimpractical for obvious reasons. On the other hand, if the patient movesslowly, only part of the detector is used. This leads to undesirableconsequences: part of the x-ray dose is wasted, discretization artifactsare enhanced, noise stability is reduced, etc.

SUMMARY OF THE INVENTION

[0007] A primary objective of the invention is to provide 3PI algorithmsfor reconstructing images of objects that have been scanned in a spiralfashion with two-dimensional detectors. For image reconstruction at anygiven voxel these algorithms require a longer section of the spiral thanthe 1PI algorithms of U.S. patent application Ser. No. 10/143,160 filedMay 10, 2002, now U.S. Pat. No. 6,574,299, which is incorporated byreference, which claims the benefit of priority to U.S. ProvisionalApplication 60/312,827 filed Aug. 16, 2001. Consequently, the newalgorithms allow to slow the patient down by about a factor of three,but still use the same size detector array.

[0008] A first preferred embodiment of the invention uses a five overallstep process for reconstructing the image of an object under a spiralscan. In a first step a current CB projection is measured. Next,families of lines are identified on a detector according to a novelalgorithm. Next, a computation of derivatives between neighboringprojections occurs and is followed by a convolution of the derivativeswith a filter along lines from the selected families of line. Next,using the filtered data, the image is updated by performing backprojection. Finally, the preceding steps are repeated for each CBprojection until an entire object has been scanned. This embodimentworks with keeping several (approximately 2-4) CB projections in memoryat a time and uses one family of lines.

[0009] For the second embodiment, different families of lines can beused in combination with keeping several CB projections in memory.

[0010] Modifications of these embodiments are possible, that will allowkeeping only one CB projection in computer memory at a time. This can bedone analogously to what was done in U.S. patent application Ser. No.10/143,160 filed May 10, 2002, which is incorporated by reference, whichclaims the benefit of priority to U.S. Provisional Application60/312,827 filed Aug. 16, 2001.

[0011] Consequently, the new algorithms allow to slow the patient downby about a factor of three, but still use the same size detector array.

[0012] Further objects and advantages of this invention will be apparentfrom the following detailed description of the presently preferredembodiments, which are illustrated schematically in the accompanyingdrawings.

BRIEF DESCRIPTION OF THE FIGURES

[0013]FIG. 1 shows a typical arrangement of a patient on a table thatmoves within a rotating gantry having an x-ray tube source and adetector array, where cone beam projection data sets are received by thex-ray detector, and an image reconstruction process takes place in acomputer with a display for the reconstructed image.

[0014]FIG. 2 shows an overview of the basic process steps of theinvention.

[0015]FIG. 3 shows stereographic projection of the spiral onto thedetector plane

[0016]FIG. 4 shows the detector plane with various projections andimportant lines

[0017]FIG. 5 shows the boundary curve

[0018]FIG. 6 shows the continuous illumination case, x/R=(0,0.25,0)

[0019]FIG. 7 shows the interrupted illumination case, x/R=(−0.5, 0, 0)

[0020]FIG. 8 shows the points where various critical events occur

[0021]FIG. 9 shows the distribution of weights inside the 5IP domain inthe case of continuous illumination

[0022]FIG. 10 shows the distribution of weights inside the 5IP domainsin the case of interrupted illumination

[0023]FIG. 11 shows the family of filtering lines parallel to L₀

[0024]FIG. 12 shows the family of filtering lines tangent to Γ_(±1)

[0025]FIG. 13 shows how to choose filtering lines depending on thelocation of {circumflex over (x)}

[0026]FIG. 14 shows the family of filtering lines tangent to Γ_(±2)

[0027]FIG. 15 shows the filtering lines and the associated constantsc_(m) in different cases: {circumflex over (x)}∈F₁ (top left panel),{circumflex over (x)}∈F₂ (top right panel), {circumflex over (x)}∈F₃ andabove Γ₁ (middle panel), {circumflex over (x)}∈G₁ (bottom left panel),{circumflex over (x)}∈G₂ (bottom right panel).

[0028]FIG. 16 shows possible locations of points s₁, s₂, s₃

[0029]FIG. 17 shows the filtering lines and the associated constantsC_(m) in different cases: {circumflex over (x)}∈F₁ (left panel),{circumflex over (x)}∈F₂ (right panel).

[0030]FIG. 18 is a three substep flow chart for finding families oflines for filtering, which corresponds to step 20 of FIG. 2.

[0031]FIG. 19 is a seven substep flow chart for preparing for filtering,which corresponds to step 30 of FIG. 2.

[0032]FIG. 20 is a seven substep flow chart for filtering, whichcorresponds to step 40 of FIG. 2.

[0033]FIG. 21 is a five substep flow chart for the first part ofback-projection, which corresponds to step 50 of FIG. 2.

[0034]FIG. 22 is a three substep flow chart for the second part ofback-projection, which corresponds to step 50 of FIG. 2.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0035] Before explaining the disclosed embodiments of the presentinvention in detail it is to be understood that the invention is notlimited in its application to the details of the particular arrangementsshown since the invention is capable of other embodiments. Also, theterminology used herein is for the purpose of description and not oflimitation.

[0036] This invention is a continuation-in-part of U.S. patentapplication Ser. No. 10/143,160 filed May 10, 2002, now U.S. Pat. No.6,574,299, which is incorporated by reference, which claims the benefitof priority to U.S. Provisional Application 60/312,827 filed Aug. 16,2001.

[0037] The invention will now be described in more detail.

Theoretical Background

[0038] First we introduce the necessary notations. Let

C:={y∈R ³ :y ₁ =R cos(s), y ₂ =R sin(s), y ₃ =s(h/2π), s∈R},  (1)

[0039] where h>0 be a spiral, and U be an open set strictly inside thespiral:

{overscore (U)}c⊂{x∈R ³ :x ₁ ² +x ₂ ² <r ²}, 0<r<R,  (2)

[0040] S² is the unit sphere in R³, and $\begin{matrix}{{{D_{f}\left( {y,\beta} \right)}:={\int_{0}^{\infty}{{f\left( {y + {\beta \quad t}} \right)}\quad {t}}}},\quad {\beta \in S^{2}},} & (3)\end{matrix}$

${{\beta \left( {s,x} \right)}:=\frac{x - {y(s)}}{{x - {y(s)}}}},{x \in U},{s \in R},$

 Π(x,ξ):={y∈R ³:(y−x)·ξ=0},  (4)

[0041] that is D_(f)(y,β) is the CB transform of f. Given(x,ξ)∈U×(R³,0), let s_(j)=s_(j)(ξ,ξ·x), j=1, 2, . . . , denote points ofintersection of the plane Π(x,ξ) with C. Also, {dot over (y)}(s):=dy/ds.For β∈S², β^(⊥) denotes the great circle {α∈S²: α·β=0}. Fix any x∈R³,where f needs to be computed. In order to compute f(x) we will use asection of the spiral of finite extent, which is to be determined later.For now it will be denoted C(x). The corresponding parametric intervalis denoted I(x). The main assumption about C(x) is the following.

[0042] Property C1. (Completeness condition) Any plane through xintersects C(x) at least at one point.

[0043] An important ingredient in the construction of the inversionformula is weight function n(s,x,α), α∈β^(⊥)(s,x). The function n can beunderstood as follows. x and α determine the plane Π(x, α), and theweight n assigned to y(s)∈Π(x,α)∩C depends on the location of x. In viewof this interpretation we assume n(s, x, α)=n(s, x, −α). The mainassumptions about n are the following.

[0044] Property W1. Normalization condition: $\begin{matrix}{{{\sum\limits_{j}{n\left( {s_{j},x,\alpha} \right)}} = {1\quad {a \cdot  \cdot {on}}\quad S^{2}}};} & (5)\end{matrix}$

[0045] Property W2. There exist finitely many C¹ functionsα_(k)(s,x)∈β^(⊥)(s,x), s∈I(x), such that n(s, x, α) is locally constantin a neighborhood of any (s, α), where s∈I(x) and α∈β^(⊥)(s,x),α∉∪_(k)α_(k)(s,x).

[0046] Denote

φ(s,x,θ):=sgn(α·{dot over (y)}(s))n(s,x,α), α=α(s,θ) ∈β^(⊥)(s,x),  (6)

[0047] Under assumptions C1, W1, and W2 the following general inversionformula is derived in the paper by A. Katsevich “A general scheme forconstructing inversion algorithms for cone beam CT,” InternationalJournal of Mathematics and Mathematical Sciences, Vol. 21, pp. 1305-1321(2003): $\begin{matrix}{{f(x)} = {{- \frac{1}{4\pi^{2}}}{\int_{I{(x)}}^{\quad}{\sum\limits_{m}{\frac{c_{m}\left( {s,x} \right)}{{x - {y(s)}}}\quad \times {\int_{0}^{2\pi}{\frac{\partial\quad}{\partial q}\quad {D_{f}\left( {{y(q)},{{\cos \quad \gamma \quad {\beta \left( {s,x} \right)}} + {\sin \quad {{\gamma\alpha}^{\bot}\left( {s,x,\theta_{m}} \right)}}}} \right)}{_{q = s}{{\frac{d\quad \gamma}{\sin \quad \gamma}{s}},}}}}}}}}} & (7)\end{matrix}$

[0048] where θ_(m)∈[0, π) are the points where θ(s, x, θ) isdiscontinuous, and c_(m)(s, x) are values of the jumps: $\begin{matrix}{{c_{m}\left( {s,x} \right)}:={\lim\limits_{ɛ->0^{+}}{\left( {{\varphi \left( {s,x,{\theta_{m} + ɛ}} \right)} - {\varphi \left( {s,x,{\theta_{m} - ɛ}} \right)}} \right).}}} & (8)\end{matrix}$

[0049] Unless n is chosen appropriately, the inversion formula is notnecessarily of the FBP type.

3PI Lines and Their Properties

[0050] Suppose that the x-ray source is fixed at y(s₀) for some s₀∈R.Since the detector array rotates together with the source, the detectorplane depends on s₀ and is denoted DP(s₀). It is assumed that DP(s₀) isparallel to the axis of the spiral and is tangent to the cylinder y₁²+y₂ ²=R² (cf. (1)) at the point opposite to the source. Thus, thedistance between y(s₀) and the detector plane is 2R (see FIG. 3).Introduce coordinates in the detector plane as follows. Let the d₁-axisbe perpendicular to the axis of the spiral, the d₂-axis be parallel toit, and the origin coincide with the projection of y(s₀). Projectstereographically the upper and lower turns of the spiral onto thedetector plane as shown in FIG. 3. This gives the following parametriccurves: $\begin{matrix}{{{d_{1}(s)} = {2R\frac{\sin \left( {s - s_{0}} \right)}{1 - {\cos \left( {s - s_{0}} \right)}}}},\quad {{d_{2}(s)} = {\frac{h}{\pi}\frac{s - s_{0}}{1 - {\cos \left( {s - s_{0}} \right)}}}},} & (9)\end{matrix}$

 ρ+2π(j−1)≦s−s₀≦2πj−ρ,j≧1, or  (10)

ρ+2πj≦s−s ₀≦2π(j+1)−ρ,j≦−1,  (11)

[0051] where ρ is determined by the radius of support of the patient:ρ=2 cos⁻¹(r/R) (cf. (2)). These curves are denoted Γ_(j), j=±1, ±2, . .. (see FIG. 4). {circumflex over (x)} denotes the projection of x.Connect an arbitrary source position y(s₀) to all points y(s) on thespiral, where

2π<s−s₀<4π or −4π<s−s₀<−2π.  (12)

[0052] This gives two surfaces, which are denoted S_(U) ^(3PI)(s₀) andS_(L) ^(3PI)(s₀). The region bounded by the two surfaces and thecylinder x₁ ²+x₂ ²=R² is denoted V_(3PI) (s₀). Let x be fixed. If s₀ issufficiently small, then S_(U) ^(3PI)(s₀) is below x. As so increases,two cases are possible. In the first one, known as continuousillumination, x enters V^(3PI()s₀) through S_(U) ^(3PI)(s₀) and leavesV^(3PI) (s₀) through S_(L) ^(3PI)(s₀). Clearly, the above procedureyields a unique 3PI interval [b₀(x), t₀(x)] and the corresponding unique3PI line L₀ ^(3PI)(x). In the second case, known as interruptedillumination, x enters and leaves V^(3PI)(s₀) several times. Moreprecisely, x intersects each of the surfaces S_(U) ^(3PI)(s₀) and S_(L)^(3PI)(s₀) exactly three times. Therefore, the above procedure now givesthree 3PI lines L_(i) ^(3PI)(x),i=1,2,3. The corresponding values of theparameter are denoted b_(i)(x), t_(i)(x), i=1,2,3. Due to the symmetryof the spiral, if x leaves (enters) V^(3PI)(s₀) when s₀=b_(i)(x), then xenters (leaves) V^(3PI)(s₀) when s₀=t_(i)(x), i=1,2,3. Consider theplane x₃=0. The boundary between the two cases is shown in FIG. 5. Wesee that the curve has no self intersections, so it divides the opendisk x₁ ²+x₂ ²<1 into two regions: X₁and X₃. In the central one, denotedX₁, there is one 3PI line for each x. In the exterior one, denoted X₃,there are three 3PI line for each x.

[0053] The following properties can be established. If there is only one3PI line for a point x, then

x∈V ^(3PI)(s₀)

s₀ ∈I _(3PI)(x):=[b₀, t₀], (13a)

2π<:t _(i)−b_(i)<4π, i=1,2,3.  (14a)

[0054] If there are three 3PI lines for a point x, then

x∈V ^(3PI)(s₀)

S₀ ∈I ^(3PI)(x):=[b ₁ ,b ₂ ]∪[b ₃ ,t ₁ ]∪[t ₂ ,t ₃],  (13b)

2π<t _(i) −b _(i)<4π, i=1,2,3.  (14b)

[0055] In terms of the detector plane equations (13a) and (13b) implythat {circumflex over (x)} is between Γ₂ and Γ⁻² if and only ifs∈I^(3PI)(x). Therefore, using the analogy with the 1PI case, we cancall this region on the detector the 3PI window. In the 1PI case, theparametric interval bounded by the endpoints of the 1PI line of x iscalled the 1PI parametric interval of x. Similarly, I^(3PI)(x) is calledthe 3PI parametric interval of x.

An Auxiliary Construction

[0056] Suppose first that continuous illumination takes place. Considerthe following two curves on the surface of the unit sphere S². Thisfirst curve consists of all unit vectors a orthogonal to L₀ ^(3PI)(x)and is denoted by A. The other curve consists of vectors $\begin{matrix}{{{\alpha (s)} = {\pm \frac{\left( {x - {y(s)}} \right) \times {\overset{.}{y}(s)}}{{\left( {x - {y(s)}} \right) \times {\overset{.}{y}(s)}}}}},\quad {s \in {I^{3{PI}}(x)}},} & (15)\end{matrix}$

[0057] and is denoted by T. It is not convenient to represent thesecurves directly on S², so they will be shown on the plane usingspherical coordinates (θ₁, θ₂) defined by

S ²

α(s)=(cos θ₁ sin θ₂,sin θ₁ sin θ₂,cos θ₂), −π≦∂

[0058] Since both α and −α define the same plane, we can restrict θ₁ toany interval of length π and “glue” the opposite boundaries using theidentification

(θ₁,θ₂)≈(θ₁+π)_(mod2π), π−θ₂).  (17)

[0059] A typical situation for x/R=(0,0.25,0) is shown in FIG. 6. It isvery convenient to think about points on S² not only as unit vectors,but also as planes. Each α∈S² corresponds to a unique plane through xwith normal vector α. This correspondence is one-to-one if vectors withopposite orientation are identified. Together A and T split S² intoseveral domains: D₁, D₂ . . . By construction, for a fixed i, the numberof points in C^(3PI)(x)∩Π(x, α) is the same for all α∈D_(i). If D_(i)contains k intersection points (IPs), it will be called a k IP domain.

[0060] Suppose now that interrupted illumination takes place. In asimilar fashion, define several curves on the surface of S₂. The firstthree curves are obtained by considering unit vectors perpendicular toeach of the three 3PI axes. They are denoted A_(k), k=1,2,3. The secondset of curves is obtained by restricting s in (15) to the intervals [b₁,b₂], [b₃,t₁], and [t₂,t₃]. These curves are denoted T₁₂, T₃₁, and T₂₃,respectively. A typical situation for x/R=(−0.5,0,0) is shown in FIG. 7.

[0061] It is clear that if x₁ and x₂ are close to each other, they willshare similar diagrams. By this we mean that by a smooth sequence oftransformations one diagram can be converted into the other, and incorresponding domains the number of IPs and their distribution over thesubintervals in I_(3PI) stays the same. An essential change is possibleonly in a neighborhood of x where a “critical event” occurs: threeboundaries intersect each other at one point on S². These points can befound numerically and the results are shown in FIG. 8. The smallestdistance from any of these points to the center of rotation isr0≈0.618R. Thus, in situations of interest in medical applications(r_(FOV)≦0.5R), only the following two cases can happen: continuousillumination as shown in FIG. 6 and interrupted illumination as shown inFIG. 7.

Construction of the Weight Function n

[0062] For any s∈I_(3PI)(x) determine the weight function n(s,x,α),α∈β^(⊥)(s,x), according to the following rule:

[0063] In 1IP domains the only IP gets weight n=1;

[0064] In 3IP domains each IP gets weight n=⅓;

[0065] In 5IP domains two IPs get weight n=0, and all the remaining IPsget weight n=⅓(see FIGS. 9 and 10).

First Reconstruction Algorithm

[0066] Once the weights are found, we have all the ingredients neededfor constructing the algorithm. First of all, for image reconstructionat x we take the 3PI spiral segment of x: C(x)=C^(3PI)(x) andI(x)=I^(3PI)(x) (cf. Section 1). As follows from (6), (7), for eachs∈I^(3PI)(x) the filtering directions are determined by finding thediscontinuities of φ(s,x,α)=sgn(α·{dot over (y)}(s))n(s,x,α). The studyof these discontinuities leads us to the following three families oflines.

[0067] The first family consists of lines parallel to the spiraltangent. This family is denoted

₀. The lines and the associated coefficients c_(m) are shown in FIG. 11.Obviously, for any {circumflex over (x)} there is only one line L∈

₀ that contains {dot over (x)}.

[0068] In FIG. 12 we see the second family of filtering lines. Itconsists of lines tangent to Γ_(±1) and is denoted

₁. Here for each {circumflex over (x)} within the 1PI window we take twolines from L₁. These lines are determined from the rule explained inFIG. 13. This figure also shows the associated coefficients c_(m).

[0069] In FIG. 14 (left panel) we see the family of filtering linestangent to Γ_(±2). This family is be denoted

₂. Now, for any {circumflex over (x)}∈F₁∪F₂∪F₄∪F₅ there might be more thone line L∈L₂ that contains {circumflex over (x)}. The unique line isdetermined from the following rule (see FIG. 14, right panel). If{circumflex over (x)}∈F₁∪F₄, the point of tangency is to the right of{circumflex over (x)}. If {circumflex over (x)}∈F₂∪F₅, the point oftangency is to the left of {circumflex over (x)}. In all cases c_(m)=⅓.

[0070]FIG. 15 summarizes the information contained in FIGS. 11-14. Itshows all possible cases where {circumflex over (x)} might be and allthe associated filtering directions and constants c_(m). In all casesthe direction of filtering as assumed to be from left to right. Itfollows from our construction (cf. FIG. 15) that θ_(m)(s,x) and c_(m)(s,x) in the general reconstruction formula (7), (8) depend on x only viaβ(s, x). Therefore, we can replace x by β(s,x) in the arguments of c_(m)and α^(⊥)and rewrite (7) in the form $\begin{matrix}{{{f(x)} = {{- \frac{1}{4\pi^{2}}}{\int_{I^{3{PI}}{(x)}}^{\quad}{\frac{1}{{x - {y(s)}}}{\overset{M{({s,\beta})}}{\sum\limits_{m = 1}}{{c_{m}\left( {s,\beta} \right)}{\Psi_{m}\left( {s,{\beta \left( {s,x} \right)}} \right)}{s}}}}}}},} & (18) \\{{\Psi_{m}\left( {s,\beta} \right)}:={\int_{0}^{2\pi}{\frac{\partial\quad}{\partial q}\quad {D_{f}\left( {{y(q)},{{\cos \quad \gamma \quad \beta} + {\sin \quad {{\gamma\alpha}^{\bot}\left( {s,\beta,\theta_{m}} \right)}}}} \right)}{_{q = s}{{\frac{d\quad \gamma}{\sin \quad \gamma}{s}},}}}}} & (19)\end{matrix}$

 x ₁ ² +x ₂ ²≦(0.618R)².  (20)

[0071] Step 31. Fix a line L from the said set of lines obtained in Step20. Note also that given a filtering line L, all x whose projectionsbelong to L and satisfy the rules mentioned above share the samefiltering line L (cf. FIG. 15). Thus, (19) becomes a convolution, (18)becomes backproj ection, and the algorithm (18), (19) is of theconvolution-based FBP type. This can be seen similarly to U.S. patentapplication Ser. No. 10/143,160 filed May 10, 2002, which isincorporated by reference. Other methods and techniques forbackprojection can be used. See additionally, for example, U.S. Pat. No.6,574,299 to Katsevich, which is incorporated by reference.

Improved Reconstruction Algorithm

[0072] Let ψ(t) be any smooth function defined on R and with theproperties ψ(0)=0, ψ′(t)>0, t∈R. Define the new family of lines

₂′ by requesting that any given line L∈

₂′ has three points of intersection with Γ_(±1)∪Γ_(±2): s₁,s₂,s₃, and thsatisfy:

s ₁ −s=ψ(s ₃ −s ₂), s+2π<s ₃ <s+4π,  (21)

s ₃ −s _(2=ψ() s ₁ −s), s−4π<s₃ <s−2π.  (22)

[0073] Recall that s is the current source position. The lines L∈

₂′ can be parameterized, for example, by s₃, 2π<|s₃|<4π. Location of s₁and s₂ depends on where s₃ is and is illustrated in FIG. 16.

[0074] Using the properties of ψ it is easy to establish that for each{circumflex over (x)}∈F_(,)∪F₂∪F₄∪F₅ there is a unique L∈

₂′ that contains {circumflex over (x)}. Also, if {circumflex over(x)}→L₂ ^(cr), then s₂,s_(3→Similarly, if {circumflex over (x)}→L) ⁻²^(cr), then s₂,s₃→−2Δ and s₁→s.

[0075] More importantly, in the case {circumflex over (x)}∈F₁∪F₂∪F₄∪F₅it is possible to reduce the number of filtering directions from two toone. The additional benefit is the improved detector usage. Thus, thetop two panels shown in FIG. 15 should be replaced by the diagrams shownin FIG. 17. The case when {circumflex over (x)} appears below L₀(i.e.,{circumflex over (x)}∈F₄∪F₅) can be obtained from FIG. 17 by symmetrywith respect to the origin in the detector plane.

[0076] The form of the inversion formula (18), (19), remains the same.The first difference is that M(s,β)=1 (instead of 2) when {circumflexover (x)}∈F₁∪F₂∪F₄∪F₅. The second difference is the direction offiltering, which is now determined from (21), (22).

General Overview of the Reconstruction Algorithms

[0077]FIG. 2 shows an overview of the basic process steps 10, 20, 30,40, 50 of the invention. The steps will now be described.

[0078] Step 10. Load the current CB (cone beam) projection into computermemory. Suppose that the mid point of the CB projections currentlystored in memory is y(s₀).

[0079] Step 20. Finding families of lines for filtering.

[0080] Using either FIG. 15 or, additionally FIG. 17 in case of theimproved algorithm, identify the required families of lines. Then,choose a discrete set of lines from each family.

[0081] Step 30. Preparing for filtering.

[0082] Parameterize points on the said lines selected in Step 20 bypolar angle. Using interpolation compute the derivative of the CB data(∂/∂q)D_(f)(y(q),β)|_(q=s) ₀ at points β on the said lines thatcorrespond to discrete values of the polar angle.

[0083] Step 40. Filtering. For each line identified in Step 20 convolvethe data for that line computed in Step 30 with filter 1/sinγ.

[0084] Step 50. Back-projection. For each reconstruction point xback-project the filtered data found in Step 40 according to equation(18). Then go to Step 10, unless there are no new CB projections toprocess or image reconstruction at all the required points x have beencompleted. Now we describe the algorithm in detail following the fivesteps 10-50 shown in FIG. 2.

[0085] Step 10. Load the current CB (cone beam) projection into computermemory. Suppose that the mid point of the CB projections currentlystored in memory is y(s₀). The detector plane corresponding to the x-raysource located at y(s₀) is denoted DP(s₀).

[0086] Step 20. Finding families of lines for filtering.

[0087] The set of lines can be selected by the following substeps 21,22, and 23.

[0088] Step 21. From the family of lines

₀ choose an equidistant set of lines that are parallel to the spiraltangent and that cover the projection of the region of interest onto thedetector plane located between Γ₂ and Γ⁻² (see FIG. 11).

[0089] Step 22. From the family of lines

₁ choose a discrete set of lines that are tangent to Γ₁ and Γ⁻¹ (seeFIG. 12). The extreme left point of tangency on Γ₁ does not have toextend beyond the projection of the region of interest onto the detectorplane. Similarly, the extreme right point of tangency on Γ⁻¹ does nothave to extend beyond the projection of the region of interest onto thedetector plane. Step 23. From the family of lines

₂ choose a discrete set of lines that are tangent to Γ₂ and Γ⁻² (seeFIG. 14, left panel). In both cases the points of tangency do not haveto extend beyond the projection of the region of interest onto thedetector plane.

[0090] In case of the improved algorithm, instead of the lines tangentto Γ₂ and Γ⁻², we choose a discrete (say, equidistant) set of values fors₃ on the curves Γ₂ and Γ⁻² and then determine the lines L∈

₂′ by solving equations (21), (22). On both curves the points s₃ do nothave to extend beyond the projection of the region of interest onto thedetector plane.

[0091] Step 30. Preparing for filtering

[0092] Step 31. Fix a line L from the said set of lines obtained in Step20.

[0093] Step 32. Parameterize points on the said line by polar angle γ inthe plane through y(s₀) and L.

[0094] Step 33. Choose a discrete set of equidistant values γ_(j) thatwill be used later for discrete filtering in Step 40.

[0095] Step 34. For each γ_(j) find the unit vector β_(j) which pointsfrom y(s₀) towards the point on L that corresponds to y_(j).

[0096] Step 35. Using the CB projection data D_(f)(y(q), Θ)) for a fewvalues of q close to so find numerically the derivative(∂/∂q)D_(f)(y(q), Θ)|_(q=s) ₀ for all θ=β_(j).

[0097] Step 36. Store the computed values of the derivative in computermemory.

[0098] Step 37. Repeat Steps 31-36 for all lines L identified in Step20. This way we will create the processed CB data corresponding to thex-ray source located at Y(s₀).

[0099] Step 40. Filtering

[0100] Step 41. Fix a line L from one of the families of lines

_(m) identified in Step 20.

[0101] Step 42. Compute FFT of the values of the said processed CB datacomputed in Step 30 along the said line.

[0102] Step 43. Compute FFT of the filter 1/sin γ

[0103] Step 44. Multiply FFT of the filter l/sin γ (the result of Steps43) and FFT of the values of the said processed CB data (the result ofSteps 42).

[0104] Step 45. Take the inverse FFT of the result of Step 44.

[0105] Step 46. Store the result of Step 45 in computer memory.

[0106] Step 47. Repeat Steps 41-46 for all lines in the said families oflines. This will give the filtered CB data Ψ_(m)(s₀,β_(j)), where mstands for the line family number from which L was selected, m=0,1,2.

[0107] By itself the filtering step is well known in the field and canbe implemented, for example, as shown and described in U.S. Pat. No.5,881,123 to Tam, which is incorporated by reference.

[0108] Step 50. Back-projection

[0109] Step 51. Fix a reconstruction point x, which represents a pointinside the patient where it is required to reconstruct the image.

[0110] Step 52. If so belongs to I^(3PI)(x), then the said filtered CBdata affects the image at x and one performs Steps 53-58. If so is notinside the interval I_(3PI)(x), then the said filtered CB data is notused for image reconstruction at x. In this case go back to Step 51 andchoose another reconstruction point.

[0111] Step 53. Find the projection {circumflex over (x)} of x onto thedetector plane DP(s₀) and the unit vector β(s₀,x), which points fromy(s₀) towards x.

[0112] Step 54. Using FIGS. 11, 13, and the right panel of FIG. 14identify the lines from the said families of lines and points on thesaid lines that are close to the said projection {circumflex over (x)}.If x is above L₂ ^(cr) or below L⁻² ^(cr), then in the case of theimproved algorithm use equations (21), (22) and FIG. 16 to find thefiltering lines close to the said projection {circumflex over (x)}. Thiswill give a few values of Ψ_(m)(s₀,β_(j)) for β_(j) close to β(s₀, x).

[0113] Step 55. With interpolation estimate the value of Ψ_(m)(s₀,β(s₀,x)) from the said values of Ψ_(m)(s₀,β_(j)) for β_(j) close to β(s₀,x).

[0114] Step 56. Compute the contribution from the said filtered CB datato the image being reconstructed at the point x by multiplying Ψ_(m)(s₀,β(s₀, x)) by −c_(m)(s,β(s₀, x))/[4π₂|x−y(s₀)|]. The appropriatebackprojection coefficient c_(m) is selected using FIGS. 11, 13, and theright panel of FIG. 14 (see also FIG. 15 for the summary). If{circumflex over (x)} is above L₂ ^(cr) or below L⁻² ^(cr), then in thecase of the improved algorithm use FIG. 17 to find the appropriatebackprojection coefficient Cm.

[0115] Step 57. Add the said contributions to the image beingreconstructed at the point x according to a pre-selected scheme (forexample, the Trapezoidal scheme) for approximate evaluation of theintegral in equation (18).

[0116] Step 58. Go to Step 51 and choose a different reconstructionpoint x. If all reconstruction points x have been processed, go to Step59.

[0117] Step 59. Go to Step 10 and load the next CB projection intocomputer memory.

[0118] The image can be displayed at all reconstruction points x forwhich the image reconstruction process has been completed (that is, allthe subsequent CB projections are not needed for reconstructing theimage at those points). Discard from the computer memory all the CBprojections that are not needed for image reconstruction at points wherethe image reconstruction process has not completed. The algorithmconcludes when the scan is finished or the image reconstruction processhas completed at all the required points.

[0119] For example, if the detector does not provide all the data whichis required for the 3PI algorithm, one can employ various techniques forestimating the missing data. In this case by the detector (respectively,cone beam projection) we mean the virtual detector (respectively,virtual cone beam projection), which includes both measured andestimated data. If one is able to estimate the missing data exactly,then the 3PI algorithm will produce exact reconstruction. In this sensewe still talk about exact reconstruction under real circumstances, whenmissing data are found approximately.

[0120] While the invention has been described, disclosed, illustratedand shown in various terms of certain embodiments or modifications whichit has presumed in practice, the scope of the invention is not intendedto be, nor should it be deemed to be, limited thereby and such othermodifications or embodiments as may be suggested by the teachings hereinare particularly reserved especially as they fall within the breadth andscope of the claims here appended.

I claim:
 1. A method of reconstructing images from data provided by atleast one detector, comprising the steps of: scanning an object in thespiral fashion with at least one detector that detects at least one conebeam projection, the cone beam projection being wider in the axialdirection than projections of four turns of the spiral that are adjacentto a current source position; and reconstructing an exact image of thescanned object in an efficient manner with a convolution based FBP(Filtered Back Projection) algorithm.
 2. The method of claim 1, whereinthe scanning step includes acquiring two-dimensional cone beam (CB)projection data of the object using the detectors.
 3. The method ofclaim 2, further comprising the step of: using the detectorssubstantially similar to those required for a 1PI algorithm.
 4. Themethod of claim 2, wherein the scanning step further includes the stepof: detecting the cone beam projection being wider in the axialdirection as compared to a cone beam projection used in a 1PI algorithm.5. The method of claim 4, wherein the scanning step further includes thestep of: detecting the cone beam projection being wider by a factor ofat least three times in the axial direction as compared to a cone beamprojection used in a 1PI algorithm.
 6. The method of claim 1, whereinthe object includes: a person.
 7. A method of computing exact imagesderived from spiral computer tomography scan with area detectors,comprising the steps of: (a) collecting cone beam (CB) projection datafrom a detector, which is wider than what is required for a 1PIalgorithm; the cone beam covering projections of four turns of thespiral that are adjacent to a current source position; (b) identifyingfamilies of lines on a plane II intersecting the cone beam projection;(c) preprocessing the CB projection data; (d) convolution-filtering saidpreprocessed CB projection data along said lines; (e) back projectingsaid filtered data to form a precursor of said image; and (f) repeatingsteps a, b, c, d, e, until an exact image of the object is completed. 8.The method of claim 7, wherein the scan includes an x-ray exposure ofthe object.
 9. The method of claim 7, wherein the steps (a)-(f) include:a 3PI algorithm.
 10. A method of computing images derived from computertomography scan with detectors, comprising the steps of: (a) collectingcone beam (CB) data from a detector during a scan of an object; (b)identifying three families of lines on a plane DP(s) intersecting thecone beam, wherein s is value of the parameter describing the scan pathand corresponding to the current source position, and the three familiesof lines include: (bi) a first family of lines parallel to {dot over(y)}(s), where {dot over (y)}(s) is the direction of the scan tangent atthe current source position; (bii) a second family of lines tangent toΓ₁ and Γ⁻¹, where Γ₁ is the projection of the scan turn defined bys<q<s+2π onto the plane DP(s); Γ⁻¹ is the projection of the scan turndefined by s−2π<q<s onto the plane DP(s); q is the parameter along thescan path which describes the point being projected; (biii) a thirdfamily of lines tangent to Γ₂ and Γ⁻², where Γ₂ is the projection of thescan turn defined by s+2π<q<s+4π onto the plane DP(s); Γ⁻² is theprojection of the scan turn defined by s−4π<q<s−2π onto the plane DP(s);(c) preprocessing and shift invariant filtering said data along saidlines of said three families; (d) back projecting said filtered data toform a precursor of said image; and (e) repeating steps a, b, c, and duntil an image of the object is completed.
 11. The method of claim 10,wherein the preprocessing includes calculation of the derivative of theCB data with respect to source position.
 12. The method of claim 10,wherein the shift invariant filtering includes convolving the saidpreprocessed data with filter 1/sin γ.
 13. The method of claim 10,wherein back projecting said filtered data from the first family oflines involves multiplying the said filtered data by the coefficientc_(m)=⅔, when the projection of x onto DP(s) is located between L₂ ^(cr)and L⁻² _(cr), where L₂ ^(cr) is the line parallel to {dot over (y)}(s)and tangent to Γ₂; L⁻² ^(cr) is the line parallel to {dot over (y)}(s)and tangent to Γ⁻².
 14. The method of claim 10, wherein back projectingsaid filtered data from lines in the first family of lines involvesmultiplying the said filtered data by the coefficient c_(m)=⅓, when theprojection of x onto DP(s) is located above L₂ ^(cr) or below L⁻² ^(cr).15. The method of claim 10, wherein back projecting said filtered datafrom a line in the second family of lines involves multiplying the saidfiltered data by the coefficient c_(m)=⅔, when the projection of x ontoDP(s) is located between Γ₁ and Γ⁻¹ and the point where the line istangent to Γ₁∪Γ⁻¹ is inside the 1PI parametric interval of x.
 16. Themethod of claim 10, wherein back projecting said filtered data from aline in the second family of lines involves multiplying the saidfiltered data by the coefficient c_(m)=−⅔, when the projection of x ontoDP(s) is located between Γ₁ and Γ⁻¹ and the point where the line istangent to Γ₁∪Γ⁻¹ is outside the 1PI parametric interval of x.
 17. Themethod of claim 10, wherein back projecting said filtered data from aline in the third family of lines involves multiplying the said filtereddata by the coefficient c_(m)=⅓.
 18. A method of computing imagesderived from computer tomography scan with detectors, comprising thesteps of: (a) collecting cone beam data from a detector during a scan ofan object; (b) identifying three families of lines on a plane DP(s)intersecting the cone beam, wherein s is value of a parameter describingthe scan path and corresponding to the current source position, and thethree families of lines include: (bi) a first family of lines parallelto {dot over (y)}(s), where {dot over (y)}(s) is the direction of thescan tangent at the current source position; (bii) a second family oflines tangent to Γ₁ and Γ⁻¹, where Γ₁ is the projection of the scan turndefined by s<q<s+2π onto the plane DP(s); Γ⁻¹ is the projection of thescan turn defined by s−2π<q<s onto the plane DP(s); q is the parameteralong the scan path which describes the point being projected; (biii) athird family of lines on the plane DP(s) that have at least three pointsof intersection s₁, s₂, s₃ with Γ_(±1) and Γ_(±2) , where Γ₂ is theprojection of the scan turn defined by s+2π<q<s+4π onto the plane DP(s);Γ⁻² is the projection of the scan turn defined by s−4π<q<s−2π onto theplane DP(s); (c) preprocessing and shift invariant filtering said dataalong said lines of said three families; (d) back projecting saidfiltered data to form a precursor of said image; and (e) repeating stepsa, b, c, and d until an image of the object is completed.
 19. The methodof claim 18, wherein the points of intersection s₁, s₂, s₃ aredetermined according to the following rules: s ₁ −s=ψ(s ₃ −s ₂) ifs+2π<s ₃ <s+4π, s ₃ −s ₂=ψ(s₁ −s) if s−4π<s ₃ <s−2π, where ψ(t) is afunction with the properties ψ(0)=0ψ′(t)>0, t∈R.
 20. The method of claim18, wherein the preprocessing includes calculation of the derivative ofthe CB data with respect to source position.
 21. The method of claim 18,wherein the shift invariant filtering includes convolving the saidpreprocessed data with filter 1/sin γ.
 22. The method of claim 18,wherein back projecting said filtered data from lines in the firstfamily of lines involves multiplying the said filtered data by thecoefficient c_(m)=⅔, when the projection of x onto DP(s) is locatedbetween L₂ ^(cr) and L⁻² ^(cr).
 23. The method of claim 18, wherein backprojecting said filtered data from a line in the second family of linesinvolves multiplying the said filtered data by the coefficient c_(m)=⅔,when the projection of x onto DP(s) is located between Γ₁ and Γ⁻¹ andthe point where the line is tangent to Γ₁∪Γ⁻¹ is inside the 1PIparametric interval of x.
 24. The method of claim 18, wherein backprojecting said filtered data from a line in the second family of linesinvolves multiplying the said filtered data by the coefficient c_(m)=−⅔,when the projection of x onto DP(s) is located between Γ₁ and Γ⁻¹ andthe point where the line is tangent to Γ₁∪_(Γ) ⁻¹ is outside the 1PIparametric interval of x.
 25. The method of claim 18, wherein backprojecting said filtered data from lines in the third family of linesinvolves multiplying the said filtered data by the coefficient c_(m)=⅔,when the projection of x onto DP(s) is located above L₂ ^(cr) or belowL⁻² ^(cr).
 26. A method of identifying a family of lines used forreconstructing images based on a scan of an object in a computertomography system, comprising the steps of (i) fixing current sourceposition s, where s is a parameter describing the scan path; (ii)picking a plane DP(s) intersecting the cone beam projection; (iii)choosing three points on the scan path, which are described using valuesof the parameter as s₁, s₂, s₃, such that: (iiia) |s−s₁<2π, (iiib)either 2π<s₂−s<4π and 2π<s₃−s<4π or −4π<s₂−s<−2π and −4π<s₃−s<−2π;(iiic) s₁−S=ψ(s₃−s₂) if 2π<s₃−s<4π or s₃−s₂=ψ(s₁−s) if −4π<s₃−s<−2π,where ψ(t) is a function with ψ(0)=0. ψ′(t)>0, t∈R; (iv) projecting thethree said points onto DP(s); and (v) drawing a line through the saidprojections.